NUWC-NPT  Technical  Report  10,840 
15  September  1997 


Evaluation  of  Small  Tail  Probabilities 
Directly  from  the  Characteristic  Function 

Albert  H.  Nuttall 

Surface  Undersea  Warfare  Directorate 


19971021  036 


Naval  Undersea  Warfare  Center  Division 
Newport,  Rhode  Island 


Approved  for  public  release;  distribution  is  unlimited. 


PREFACE 


The  work  described  in  this  report  was  sponsored 
by  the  Independent  Research  Program  of  the  Naval 
Undersea  Warfare  Center  Division,  Newport,  RI , 
under  Project  No.  B100077,  "Near-Optimum  Detection 
of  Random  Signals  with  Unknown  Locations, 

Structure,  Extent,  and  Strengths,"  principal 
investigator  Albert  H.  Nuttall  (Code  311).  The 
independent  research  program  is  funded  by  the 
Office  of  Naval  Research;  the  Naval  Undersea 
Warfare  Center  Division  Newport  program  manager  is 
Stuart  C.  Dickinson  (Code  102). 


The  technical  reviewer  for  this  report  was 
Paul  M.  Baggenstoss  (Code  2121). 


Reviewed  and  Approved:  15  September  1997 


Director,  Surface  Undersea  Warfare 


REPORT  DOCUMENTATION  PAGE 

Form  Approved 

0MB  No.  0704-0188 

Public  reporting  burden  for  this  colleaion  of  information  is  estimated  to  average  l  hour  per  response,  including  the  time  for  reviewing  instruaions.  searching  existing  data  sources, 
gathering  and  maintaining  the  data  needed,  and  completing  and  reviewing  the  colleaion  of  information.  Send  comments  reoarding  this  burden  estimate  or  any  other  aspea  of  this 
colleaion  of  information,  including  suggestions  for  r^ucing  this  burden,  to  Washington  Headauaaers  Services,  Direaorate  for  information  Operations  and  Reports,  1215  Jefferson 
Davis  Highway,  Suite  1204,  Arlington.  VA  22202-4302,  and  to  the  Office  of  Management  and  Budget.  Paperwork  Reduaion  Projea  (0704-0 1B8).  Washington,  DC  20503. 

1.  AGENCY  USE  ONLY  (Leave  blank)  2.  REPORT  DATE  3.  REPORT  TYPE  AND  DATES  COVERED 

15  September  19S7  Proaress 

4.  TITLE  AND  SUBTITLE 

Evaluation  of  Small  Tail  Probabilities 

Directly  from  the  Characteristic  Function 

5.  FUNDING  NUMBERS 

PE  0601152N 

6.  AUTHOR(S) 

Albert  H.  Nuttall 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Naval  Undersea  Warfare  Center  Division 

1176  Howell  Street 

Newport,  Rhode  Island  02841-1708 

8.  PERFORMING  ORGANIZATION 

REPORT  NUMBER 

NUWC-NPT  TR  10,840 

9.  SPONSORING /MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES)  ' 

Office  of  Naval  Research 

800  North  Quincy  Street,  BCT  1 

Arlington,  VA  22217-5000 

10.  SPONSORING /MONITORING 

AGENCY  REPORT  NUMBER 

11.  SUPPLEMENTARY  NOTES 

12a.  DISTRIBUTION /AVAILABILITY  STATEMENT 

Approvecd  for  public  release; 
distribution  is  unlimited. 

12b.  DISTRIBUTION  CODE 

13.  ABSTRACT  (Maximum  200  words) 

.  An  ef f i ci ent /  fast,  and  accurate  Fourier  transform  technique  for 
obtaining  small  tail  probabilities  for  both  the  probability  density 
function  and  the  exceedance  distribution  function,  directly  from  the 
characteristic  function,  is  derived  and  demonstrated  numerically  for 
several  examples.  The  method  is  especially  useful  when  analytic  or 
asymptotic  expressions  for  the  probabilities  are  unavailable  or 
unknown . 

By  choosing  the  shift  parameter  r  close  to  the  highest 
singularity  of  the  characteristic  function  in  the  complex  plane, 
very  small  values  of  the  tail  probabilities  of  the  density  function 
and  exceedance  distribution  function  can  be  realized.  The  cost  in 
this  approach  is  that  the  sampling  increment  must  then  be  taken 
small,  in  order  to  avoid  aliasing.  Finer  sampling  necessitates  more 
computer  time  and  effort,  but  it  does  not  require  more  storage; 

14.  SUBJECT  TERMS 

Probabilities  Characteristic  Function 
Probability  Density  Exceedance  Distribution 
Fast  Fourier  Transform  Aliasina 

15.  NUMBER  OF  PAGES 

54 

16.  PRICE  CODE 

17.  SECURITY  CLASSIFICATION 

OF  REPORT 

Unclassified 

18.  SECURITY  CLASSIFICATION  ' 
OF  THIS  PAGE 

Unclassified _ 

19.  SECURITY  CLASSIFICATION 

OF  ABSTRACT 

20.  LIMITATION  OF  ABSTRACT 

- — - - - 

MSN  7540-01-280-5500  Standard  Form  298  (Rev,  2-89) 


Prpscrib^  bv  AIMS!  Std  Z39-18 


UNCLASSIFIED _ 

SECURITY  CLASSIFICATION 

OF  THIS  PAGE _ _ _ 

rather,  prealiasing  can  be  advantageously  employed  to  keep  the 
fast  Fourier  transform  size  at  reasonable  values.  The  fast 
Fourier  transform  size  has  no  effect  upon  the  errors  caused  by 
aliasing  and  truncation;  rather,  the  size  merely  controls  the 
spacing  at  which  the  probability  density  function  and 
exceedance  distribution  function  outputs  are  calculated. 

Tail  probabilities  in  the  E-50  range  are  readily  available 
with  a  computer  limited  to  15  significant  decimal  digits. 


UNCLASSIFIED _ 

SECURITY  CLASSIFICATION 
OF  THIS  PAGE 


TABLE  OF  CONTENTS 


Page 

LIST  OF  ILLUSTRATIONS . ii 

LIST  OF  ABBREVIATIONS  AND  SYMBOLS . iii 

INTRODUCTION  .  1 

DEVELOPMENT  OF  TECHNIQUE  FOR  PDF . 3 

Movement  of  Contour  for  PDF . 3 

Sampling  and  Aliasing  for  PDF . 6 

FFT  Considerations  for  PDF . 8 

Plotting  Procedure  .  10 

DEVELOPMENT  OF  TECHNIQUE  FOR  EDF . 13 

Movement  of  Contour  for  EDF  . . 13 

Sampling  and  Aliasing  for  EDF . 15 

FFT  Considerations  for  EDF . 17 

EXAMPLES . 19 

Chi-Square  PDF . 19 

Gaussian  PDF . 21 

Branch  Point  PDF . 23 

Essential  Singularity  PDF  .  25 

Chi-Square  EDF . 28 

Gaussian  EDF . 30 

Branch  Point  EDF . 30 

Essential  Singularity  EDF  . 33 

SUMMARY.  . . 35 

REFERENCES . 36 

APPENDIX  A  -  BASIC  PROGRAM  FOR  PDF . •  •  •  -A-l 

APPENDIX  B  -  BASIC  PROGRAM  FOR  EDF . B-1 


1 


LIST  OF  ILLUSTRATIONS 


Figure 

1  Chi-Square  p(u)  via  Standard  FFT  Technique 

2  Chi-Square  p(u)  . 

3  Chi-Square  p(u)  . 

4  Gaussian  p(u)  . 

5  Gaussian  p(u)  . 

6  Branch  Point  p(u)  . 

7  Branch  Point  p(u)  . 

8  Essential  Singularity  p(u)  . 

9  Essential  Singularity  p(u)  . 

10  Chi-Square  E(u)  . 

11  Chi-Square  E(u) . . 

12  Gaussian  E(u)  . 

13  Gaussian  E(u) . . 

14  Branch  Point  E(u)  . 

15  Branch  Point  E(u) . . 

16  Essential  Singularity  E(u)  . 

17  Essential  Singularity  E(u)  . 


Page 

20 

22 

22 

24 

24 

26 

26 

27 

27 

29 

29 

31 

31 

32 
32 
34 
34 


11 


LIST  OF  ABBREVIATIONS  AND  SYMBOLS 


C 

CF 

EDF 

E(u) 

E(u) 

E(u) 

f(^) 

f_(^) 

FFT 

Im 

N 

PDF 

Prob 

p{u) 

£(u) 

p(u) 

r 

Re 


Real  part  of  location  of  a  singularity 
Imaginary  part  of  location  of  a  singularity 
Minimum  imaginary  part  of  singularity  locations 
Additive  real  constant 
Contour  of  integration,  equation  (17) 

Characteristic  function,  equation  (1) 

Exceedance  distribution  function,  equation  (17) 
Exceedance  distribution  function  of  x,  equation  (17) 
Auxiliary  function,  equation  (19) 

Periodic  function,  equation  (21) 

Characteristic  function  of  x,  equation  (1) 
Characteristic  function  of  positive-u  part  of  p(u) 
Characteristic  function  of  negative-u  part  of  p(u) 
Finite  sequence  of  samples,  equation  (12) 

Fast  Fourier  transform,  equation  (13) 

Imaginary  part,  equation  (2) 

Size  of  FFT,  equation  (13) 

Integer  in  the  noise  region,  equation  (25) 
Probability  density  function,  equation  (1) 
Probability 

Probability  density  function  of  x,  equation  (1) 
Auxiliary  function,  equation  (6) 

Periodic  function,  equation  (7) 

Real  shift  in  path  of  I  integration,  equation  (4) 
Real  part,  equations  (9),  (10),  (21) 


iii 


LIST  OF  ABBREVIATIONS  AND  SYMBOLS  (Cont'd) 


General  sampling  point,  equation  (15) 

X  Random  variable 

Increment  in  u,  equation  (14) 

Sampling  increment,  equations  (7),  (21) 

Weighting  sequence,  equations  (10),  (21) 

bold  Random  variable 


IV 


EVALUATION  OF  SMALL  TAIL  PROBABILITIES 
DIRECTLY  FROM  THE  CHARACTERISTIC  FUNCTION 

INTRODUCTION 

Given  the  characteristic  function  (CF)  f(^)  of  a  random 
variable  x,  the  corresponding  probability  density  function  (PDF) 
p(u)  can  be  found  analytically  by  the  Fourier  transform 

00 

p(u)  »  J  d^  exp(-iu^)  f(5)  .  (1) 

_00 

The  CF  f(^)  exists  for  all  real  ^  because  the  area  under 
nonnegative  PDF  p(u)  is  unity.  However,  when  integral  (1)  cannot 
be  accomplished  in  closed  form,  it  is  often  necessary  to  resort 
to  a  numerical  technique,  namely,  the  fast  Fourier  transform 
(FFT)  to  obtain  approximate  values  of  PDF  p(u). 

Similarly,  the  exceedance  distribution  function  (EDF), 

00 

E(u)  =  Prob(x  >  u)  =  J  dt  p(t),  corresponding  to  CF  f(^),  is 

u 

available  by  the  Fourier  transform  (reference  1,  equation  (4.14) 
or  reference  2,  page  3) 

00 

E(u)  =7  +  ^1^  Im^f(^)  exp(-iu^)j  .  (2) 

0 

If  this  latter  integral  (2)  cannot  be  done  in  closed  form,  or  if 
the  indefinite  integral  of  PDF  p(u)  is  unavailable,  an  efficient 
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numerical  technique  employing  the  FFT  can  again  be  utilized 
(references  2  through  4). 

However,  there  is  an  inherent  limitation  in  the  numerical 
approximations  to  integrals  (1)  and  (2)  in  their  current  forms, 
as  integrals  along  the  real  ^-axis.  Both  require  sampling  and 
summing  integrands  with  values  in  the  neighborhood  of  unity  and 
relying  on  the  cancellation  of  complex  plus-and-minus  values 
(contributed  by  the  oscillating  exponential  exp(-iu^))  to  reach 
very  small  values  on  the  tails  of  the  PDF  p(u)  or  EDF  E(u)  for 
large  u.  If  the  computer  being  used  is  limited  to  D  significant 
decimal  digits,  then  the  small  tail  probabilities  calculated  for 
p(u)  or  E(u),  which  are  less  than  approximately  10  are  useless 
because  these  values  are  buried  in  the  inherent  round-off  noise. 

In  this  report,  an  efficient  technique  for  obtaining  small 
tail  probabilities,  significantly  less  than  10  is  presented 
for  both  the  PDF  p(u)  and  the  EDF  E(u).  The  basic  idea  is  to 
move  the  contour  of  integration  into  the  complex  ^-plane; 
however,  care  must  then  be  taken  to  control  the  aliasing  that 
always  accompanies  equispaced  sampling,  which  in  turn,  is 
required  for  efficient  use  of  an  FFT.  Examples  of  the  proposed 
technique  are  presented  for  both  the  PDF  and  the  EDF.  Values  of 
small  tail  probabilities  in  the  10  range  are  readily  achieved 
using  a  computer  with  D  =  15  significant  decimal  digits. 
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DEVELOPMENT  OF  TECHNIQUE  FOR  PDF 

MOVEMENT  OF  CONTOUR  FOR  PDF 

The  CF  f(^)  follows  from  the  PDF  p(u)  according  to  the 
inverse  Fourier  transform 

GO 

f(^)  =  J  du  exp(i^u)  p(u) 

^00 

0  C® 

=  J  du  exp(i^u)  p(u)  +  J  du  exp(i^u)  p(u)  =  f_(^)  +  f^(^)  .  (3) 
—  00  0 

The  component  f_(^),  corresponding  to  u  <  0 ,  is  analytic  for 

<  0  because  exp(-5^u)  is  bounded  for  u  <  0.  The  decay  of  p(u) 
as  u  ^  is  arbitrary. 

It  is  presumed  that  the  PDF  p(u)  has  already  been  shifted  on 
the  u-axis  so  that  it  is  virtually  zero  for  u  <  0,  yet  is  packed 
as  closely  as  possible  to  the  origin  at  abscissa  values  near 
u  =  0+.  This  condition  can  be  accomplished  by  adding  a  real 
constant  c  to  the  random  variable  under  investigation.  Addition 
of  a  constant  c  corresponds  to  multiplication  of  the 
corresponding  CF  by  the  factor  exp(ic^);  this  factor  is  presumed 
to  be  present  in  the  given  CF  f(^).  The  constant  c  can  be 
positive  or  negative,  depending  on  the  desired  direction  of  shift 
of  a  given  PDF,  so  as  to  make  p(u)  essentially  nonzero  only  in 
the  positive  neighborhood  of  u  =  0. 
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This  section  focuses  on  the  positive  tail  of  PDF  p(u).  It  is 
presumed  that  p(u)  decays  according  to  u'^  exp(-3u)  as  u  +<*>, 

3  >  0.  Then,  the  locations  “  ®n  ^  ^n  singularities 

of  f, (^)  must  satisfy  b„  >  0  for  all  n.  Let  b,  =  min{b  }, 
without  loss  of  generality.  Then,  p(u)  exp(ru)  decays 
exponentially  as  u  +®,  provided  that  real  constant  r  <  b^ . 
Accordingly,  the  Fourier  relation  in  equation  (1)  can  be  written 
as 

»-ir 

P(u)  =  2n  J  exp(-iu^)  f(^)  ,  0  <  r  <  b^  .  (4) 

-®-ir 

This  path  of  integration  passes  through  the  horizontal  strip  of 
analyticity  between  the  real  ^-axis  and  the  highest  singularity 
of  f(^)  in  the  lower-half  ^-plane. 

For  example,  PDF  p(u)  =  1.5  exp(-u)  -  exp(-2u),  for  u  >  0 
has  CF  f(^)  =  1.5/(l-i^)  -  l/(2-i^),  with  poles  at  ^  =  -i  and 

K  =  -i2.  Thus,  a^  =  0,  b^  =  1 ,  a2  =  0,  b2  =  2,  and  r  <  1.  As  a 

2 

second  example,  p(u)  =  u  exp(-u),  for  u  >  0  has  f(K)  =  l/(l-i^) 
with  a^^  =  0,  b^  =  1,  for  which  r  <  1.  Finally,  Gaussian  PDF 
p(u)  =  (2n)  ^  exp[ - ( u-c ) ^/2 ]  has  CF  f(^)  =  exp(-^^/2  +  ic^), 
which  has  no  singularities  anywhere  in  the  finite  ^-plane.  Then, 
r  can  be  chosen  arbitrarily,  positive  or  negative.  (Positive 
constant  c  is  taken  as  small  as  possible,  subject  to  PDF  value 
p(0)  being  sufficiently  small,  as  discussed  earlier.) 
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(5) 


Upon  making  the  substitution  x  =  ^  +  ir  in  equation  (4), 
it  follows  that 

p(u)  =  exp(-ru)  2(u)  , 
where  auxiliary  function 

GO 

2(u)  s  ^  J  dx  exp(-iux)  f(x  -  ir)  .  (6) 

_G0 

The  real  function  £(u)  is  nonnegative  for  all  u,  but  it  may  not 
have  unit  area.  The  relations  in  equations  (5)  and  ( 6 )  are 
exact,  even  if  the  path  of  integration  is  very  close  to  the 
highest  singularity  of  f(x  -  ir)  in  the  lower-half  x-plane  at 

x^  =  aj^  -  i  ( b^  -  r )  . 

By  writing  equation  (5)  in  the  form  ^(u)  =  p(u)  exp(ru),  it 
is  seen  that  the  closer  r  is  taken  to  b^^,  the  less  decay  there  is 
in  £(u)  for  positive  u.  This  property  makes  £(u)  a  wider 
function  of  u  and  brings  some  of  the  originally  very  small  values 
of  p(u)  up  to  moderate  levels;  of  course,  the  exponential  decay 
of  £(u)  will  eventually  take  over,  but  it  will  be  shifted  to 
larger  u  values  before  becoming  dominant.  Thus,  larger  values  of 
r  (nearer  b^^ )  enable  investigation  of  deeper  tails  of  p(u)  than 
possible  when  r  =  0,  namely,  using  the  real  ^-axis  in  equation 
(4)  . 

Because  PDF  p(u)  =  0  for  u  <  0,  the  same  is  true  for  £(u), 
according  to  equation  (5).  Also,  the  essential  nonzero  region  of 
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£(u)  is  in  the  positive  neighborhood  of  u  =  0. 

(If  the  negative  tail  of  random  variable  x  is  of  interest, 
consider  the  variable  y  =  T  -  x,  for  which  the  PDF  is 
Py,(u)  =  Pjj(T  -  u).  Real  constant  T  is  taken  just  large  enough 
that  y  is  essentially  nonzero  only  in  the  positive  neighborhood 
of  zero.  Then,  the  CF  of  interest  is  fy(^)  =  exp(i5T)  f^(-5) 
instead  of  f^(^).) 


SAMPLING  AND  ALIASING  FOR  PDF 

Fourier  relation  (6)  for  widened  function  £(u)  is  evaluated 
approximately  by  sampling  at  increment  for  all  x.  (Truncation 
errors  are  discussed  below.)  The  result  is 

A 

p(u)  =  >  !  exp(-iumA  )  f(mA  -  ir)  s  p(u)  for  all  u  .  (7) 

/a  Jl  XX 

in=— 

The  latter  real  function,  p(u),  is  periodic  in  u,  with  period 
2ii/A^.  In  fact,  it  is  the  aliased  version  of  £(u): 

00 

p(u)  ==  >  !  £ru  -  m  for  all  u  .  (8) 

For  p(u)  to  be  a  good  approximation  to  £(u)  in  the  positive 
neighborhood  of  u  =  0,  the  aliasing  lobes  of  p(u)  must  be 
sufficiently  separated  so  that  the  fundamental  period,  interval 
(0,2ii/A^),  is  only  slightly  contaminated  by  the  undesired  lobes 
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contributed  when  ni  =  +1,  +2,....  (Recall  that  jd(u)  is 
essentially  nonzero  only  in  the  positive  neighborhood  of  u  =  0.) 
This  requirement  can  be  checked  by  looking  at  one  period  of  p(u), 
namely,  0  <  u  <  to  determine  if  the  skirts  near  u  =  0+  and 

u  =  are  sufficiently  small.  If  not,  sampling  increment 

A  must  be  decreased.  Then,  p(u)  will  be  a  good  approximation 

X 

to  £(u)  in  the  interval  (0,2ii/A^). 

It  should  be  recalled  that  the  closer  r  is  taken  to  b^ ,  the 
less  decay  there  is  in  £(u)  (see  equation  (5));  that  is,  £(u)  is 
wider  in  u,  leading  to  more  severe  aliasing  problems  in  equation 
(8),  unless  A^  is  additionally  decreased,  separating  the  lobes  of 
p(u)  farther  apart.  This  additional  decrease  is  an  unavoidable 
consequence  of  investigating  deeper  tails  of  p(u)  by  relations 
( 5 )  and  ( 6 ) . 

The  proximity  of  the  path  of  integration  in  equation  (6) 
to  the  nearest  singularity  of  f(x  -  ir)  does  not  require  an 
additional  decrease  in  sampling  increment  A^;  it  is  all  accounted 
for  in  the  widening  effect  of  exp(ru)  upon  p(u)  in  the  function 
£(u)  =  exp(ru)  p(u)  and  the  attendant  more  stringent  alias- 
suppression  requirement  in  equation  (8). 

A  shortcut  for  the  evaluation  of  equation  (7)  is  available 
by  taking  advantage  of  the  conjugate  symmetry  of  f(x  -  ir). 

In  particular,  from  equation  (3),  it  follows  that  f(-  x  -  ir)  = 
f*(x  -  ir).  Then,  equation  (6)  can  be  expressed  as 
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£(u)  =  ^  Re  J  dx  exp(-iux)  f(x  -  ir)  ,  (9) 

0 

and  equation  (7)  can  be  expressed  as 
A  ® 

p(u)  2  ^  Re  >  i  e  exp(-iumA  )  f(niA  -  ir)  =  p(u)  .  (10) 

TL  _  A  lu  X  X 

m=U 

This  function  is  the  same  function  p(u)  encountered  in  equation 
(7).  The  sequence  {Sj^}  is  h  for  m  =  0,  and  1  otherwise. 

FFT  CONSIDERATIONS  FOR  PDF 

Periodic  function  p(u),  defined  by  equation  (7),  needs  to  be 
evaluated  over  one  period  only.  In  particular,  because  the 
function  of  interest,  £(u),  is  essentially  nonzero  only  in  the 
positive  neighborhood  of  u  =  0,  attention  is  confined  to  values 

A  ” 

I  ^  \  '  exp(-i2iimn/N)  f(mA  -  ir)  for  0  <  n  <  N-1.  (11) 

m=— ® 

However,  because  the  exponential  factor  exp(-i2Timn/N)  treats  all 
complex  samples  f({m  +  kN}A^  -  ir)  equally,  regardless  of  the 
value  of  integer ^k,  these  samples  can  be  collapsed  (or 
prealiased)  into  a  set  of  just  N  values  according  to 

00 

5  s  )  I  ff(m  +  kN)A  -  irl  for  0  <  m  <  N-1  .  (12) 

k=-®  V  X  ; 
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This  process  accounts  for  all  the  samples  of  f  encountered  in 
equation  (11).  FFT  size  N  can  be  much  smaller  than  the  number  of 
terms  required  in  equation  (12)  to  control  and  minimize 
truncation  errors.  Then^  equation  (11)  can  be  written  exactly  as 
finite  sum 


p(|^  “  2n  ^  exp(-i2iimn/N)  for  0  <  n  <  N-1 


Because  the  N  infinite  sums  required  by  equation  (12)  can 
never  be  conducted  in  practice,  it  is  necessary  to  truncate  them 
where  |f(x  -  ir)|  is  sufficiently  small.  The  effect  of  the 
resultant  truncation  errors  can  be  observed  in  the  u  domain  as  a 
low-level  oscillatory  behavior  when  equation  (13)  is  plotted.  I 
this  level  is  excessive,  a  smaller  tolerance  must  be  set  on 
truncating  equation  (12),  and  additional  values  for  larger  |k| 
must  be  taken  into  account. 


Operation  (13)  can  be  accurately  and  efficiently  realized  as 
an  N-point  FFT  if  N  is  taken  as  a  power  of  2.  Observe  that  FFT 
size  N  does  not  affect  the  truncation  error  or  the  aliasing 
error;  rather,  N  merely  sets  the  spacing  in  variable  u  in 
equation  (13),  namely. 


(14) 


at  which  output  p(u)  is  obtained. 
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PLOTTING  PROCEDURE 


It  is  recommended  that  output  (13)  be  plotted  over  the  full 
sweep  of  n  values  and  observed  for  its  near-origin  behavior, 
aliasing,  and  truncation  errors.  Then,  additive  constant  c,  or 
sampling  increment  or  the  truncation  limits  in  equation  (12), 

or  any  combination  of  the  three  can  be  modified,  and  the  entire 
procedure  repeated  until  all  errors  are  acceptable.  Also, 
regions  of  u  where  p(u)  is  below  the  inherent  round-off  error  of 
the  computer  will  appear  with  a  white— noise  strip  of  low-level 
values;  this  white-noise  strip  gives  a  measure  of  the  relative 
errors  of  the  values  obtained  for  p(u)  in  the  entire  region 
0  <  u  <  2ii/Aj^. 

If  the  following  points  in  the  fundamental  interval  are 
defined  according  to 

u  =  n  A  =  §  for  0  <  n  <  N-1  ,  (15) 

n  u  A  N  -  “ 

X 

then,  the  desired  values  of  pdf  p(u)  are  available  from  equations 
(5),  (7),  and  (13)  according  to 

p(u^)  =  exp(-rUj^)  2(u^)  =  exp(-ru^)  p(u^) 

A  N-1 

=  exp(— ru  )  5  !  exp ( —  i2 nmn/N )  f  for  0  <  n  <  N— 1  .  (16) 

^  n  2ii  i  m 

m=0 

A  plot  of  this  approximation  to  p(u),  in  conjunction  with  the 
relative  error  information  obtained  from  the  earlier  plot  of  p  in 
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equation  (13),  reveals  the  levels  of  PDF  p(u)  and  the  relative 
errors  associated  with  different  regions  of  u.  If  the  lowest 
levels  reached  by  equation  (16)  are  not  sufficiently  small, 
parameter  r  in  equations  (4)  and  (12)  must  be  increased  and  the 
plots  of  equations  (13)  and  (16)  must  be  repeated.  As  noted 
earlier,  doing  so  will  likely  entail  a  decrease  in  sampling 
increment  to  control  aliasing  of  a  wider  p(u) .  A  short 
trial-and-error  procedure  may  be  required. 


11/(12  blank) 


DEVELOPMENT  OF  TECHNIQUE  FOR  EDF 


MOVEMENT  OF  CONTOUR  FOR  EDF 

The  Fourier  transform  relating  the  EDF  E(u)  =  Prob(x  >  u)  to 
the  CF  f(^)  can  be  written  in  the  form 

00 

E(u)  =  J  dt  p(t)  =  J  exp(-iu^)  ,  (17) 

u  C 

where  contour  C  is  the  real  ^-axis,  except  that  it  passes  below 
the  singularity  at  ^  =  0,  The  comments  regarding  the 
corresponding  PDF  p(u)  in  equation  (3)  and  sequel  are  directly 
relevant  to  equation  (17).  In  particular,  as  in  equation  (4), 
the  contour  is  moved  to  distance  r  below  the  real  ^-axis,  where 
real  constant  r  is  less  than  the  distance  to  the  nearest 
singularity  of  f(^)  in  the  lower— half  plane  to  obtain 

00-  i  r 

E(u)  =  J  d^  I '  exp(-iu^)  =  exp(-ru)  E(u)  ,  (18) 

_oo_i  r 

where  auxiliary  function 


S'"’  ■  ik  J 


f(x  -  ir) 


exp( -iux ) 


Here,  the  substitution  x  =  ^  +  ir  was  used.  Relation  (19)  is 
still  a  Fourier  transform,  although  with  a  modified  integrand. 
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The  real  function  E(u)  is  nonnegative  for  all  u.  Relations 
(18)  and  (19)  are  exact,  even  if  the  path  of  integration  is  very 
close  to  the  highest  singularity  of  f(x  -  ir)  in  the  lower-half 
x-plane  at  Xj^  =  a^  -  i(b^  -  r),  >  0. 

By  writing  equation  (18)  in  the  form  E(u)  =  E(u)  exp(ru),  it 
is  seen  that  the  closer  r  is  taken  to  b^,  the  less  decay  there  is 
in  E(u)  for  positive  u.  For  u  >  0,  E(u)  becomes  a  wider  function 
of  u  and  brings  some  of  the  originally  very  small  values  of  E(u) 
up  to  moderate  levels;  of  course,  the  exponential  decay  of  E(u) 
as  u  -»  +®  eventually  takes  over,  but  it  is  shifted  to  larger  u 
values  before  becoming  dominant.  Thus,  larger  values  of  r 
(nearer  b^)  enable  investigation  of  deeper  tails  of  E(u)  than 
possible  when  r  =  0,  namely,  using  the  real  ^-axis  in  equation 
(17)  . 


The  dominant  asymptotic  behaviors  of  functions  E(u)  and  E(u) 
are  as  follows: 


'  1  as  u 

'  exp(  ru)  as  u  ^  -®' 

exp(-b^u)  as  u  +®J 

>  ,  E  ( U )  «  < 

> 

[^exp  (  ru-b^u )  as  u 

(20) 


By  taking  0  <  r  <  b^ ,  E(u)  decays  to  zero  for  both  u  ±®. 

Taking  r  close  to  b^  keeps  E(u)  at  a  high  level  (near  E(0))  for 
an  extended  range  of  positive  u.  It  also  gives  maximum  decay  to 
the  left  tail  of  E(u),  which  will  be  seen  to  be  advantageous  in 
terms  of  alias  control.  If  E(u)  can  be  accurately  calculated  for 
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large  positive  u,  then  relation  (18)  will  result  in  accurate 
evaluation  of  E(u)  at  very  low  levels  of  probability. 

It  is  important  to  observe  that,  unlike  the  original  EDF 
E(u),  the  modified  EDF  E(u)  decays  to  zero  as  u  +“>;  the 
original  EDF  E(u)  always  approaches  1  as  u  This  feature  of 

E(u)  allows  the  integral  of  equation  (19)  to  be  approximated  by  a 
sampling  procedure  and  yet  have  controllable  aliasing  lobes. 

SAMPLING  AND  ALIASING  FOR  EDF 

Fourier  relation  (19)  for  modified  EDF  E(u)  is  evaluated 
approximately  by  sampling  at  increment  for  all  x.  The  result 
is 

00 

E(u)  =  ^  Re  J  dx  +"ix^  exp(-iux)  = 

0 

A  <o  f  ( mA  -  ir ) 

'  ^  Re  C  r  I  ima  exp(-iuta4^)  .  E(u)  toe  all  u  .  (21) 

in=0  X 

The  latter  real  function,  E(u),  is  periodic  in  u,  with  period 
2n/A^;  in  fact,  it  is  the  aliased  version  of  E(u): 

00 

E(u)  =  >  !  sfu  -  m  for  all  u  .  (22) 

m=-“  ^  ^x"’ 

For  E(u)  to  be  a  good  approximation  to  E(u)  in  the  positive 
neighborhood  of  u  =  0,  the  aliasing  lobes  of  E(u)  must  be 
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sufficiently  separated  so  that  the  fundamental  period,  interval 
(0,2n/A^),  is  only  slightly  contaminated  by  the  undesired  lobes 
contributed  when  m  =  +l,+2,....  (Recall  that  E(u)  is  essentially 
nonzero  only  in  the  neighborhood  of  u  =  0;  see  equation  (20).) 
This  aliasing  requirement  can  be  checked  by  looking  at  one  period 
of  E(u),  namely,  0  <  u  <  2ii/A^,  to  see  if  the  skirts  near  u  =  0+ 
and  u  =  2ii/A^-  are  sufficiently  small.  If  they  are  not,  sampling 
increment  A^  must  be  decreased.  Then,  E(u)  can  be  made  a  good 
approximation  to  E(u)  in  an  interval  of  length  (0,2ii/A^)  located 
about  u  =  0 . 

From  equations  (18)  and  (21),  interest  is  centered  on 

E(u)  =  exp(-ru)  E(u)  =  exp(-ru)  E(u)  (23) 

for  u  in  the  fundamental  period  of  length  2ii/A^  in  the 
neighborhood  of  zero.  Accurate  results  for  EDF  E(u)  are  obtained 
when  E(u)  has  experienced  insignificant  aliasing,  which  is 
accomplished  by  choosing  A^  small  enough.  Because  E(u)  for  u  >  0 
has  purposely  been  stretched  out  over  a  wider  u  interval,  by 
choice  of  r,  A  must  be  made  quite  small.  Making  A  quite  small, 
however,  is  an  inherent  requirement  for  investigating  very  low 
levels  of  E(u)  because  these  low  levels  are  not  realized  until  a 
considerable  interval  of  u  has  been  covered.  In  other  words,  the 
more  stringent  requirement  on  A^  is  not  due  to  the  particular 
procedure,  but  rather,  due  to  the  desire  to  investigate  a  wide r 
range  of  positive  u  values,  namely,  those  including  very  small 
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E(u).  That  procedure  naturally  requires  a  larger  spacing  of  the 
aliasing  lobes  of  E(u)  to  clear  out  a  wider  interval  of  u  space 
so  that  the  ever-present  aliasing  does  not  become  intolerable. 


FFT  CONSIDERATIONS  FOR  EDF 

Because  E(u)  has  period  2ii/A^,  its  evaluation  can  be  limited 
to  one  period.  In  particular,  consider 


»  f  ( mA  -  i  r  ) 

Re  5  '  e  - - — : — r —  exp ( -i2nmn/N)  for  0  <  n  <  N-1 . 

^ — i,  m  r  +  imA  ^ 

m=0  X 


Presuming  that  sampling  increment  A^  is  taken  small  enough  to 
effectively  eliminate  aliasing,  a  plot  of  this  complete  period  of 
E(u)  reveals  an  interior  region  where  the  computer  round-off 
noise  is  dominant. 


If  an  integer  in  this  noise  region  is  defined  as  N^,  then  the 
values  of  equation  (24)  in  the  range  (0,Nj)  can  be  used  directly 
to  get  estimates  of  the  desired  EDF  according  to 


(2n 

(  2n 

n') 

„  (2n 

n't 

(  ^  2n 

s  f2ii 

n'l 

[a 

X 

nJ 

=  exp[-r  Y- 

X 

nJ 

5  — 

X 

nJ 

=  exp[-r  j- 
X 

nJ 

X 

nJ 

for  0  <  n  <  Nj.  However,  the  values  in  the  region  <  n  must  be 
scaled  differently  because  they  should  represent  the  negative 
region  of  E(u)  just  below  the  origin,  that  is,  u  <  0.  The 
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correct  scaling  in  this  region  is 


exp(r  |S-  r  E  (fl  a]  for  «,  <  n  <  N  .  (26, 

X  X  X 

These  scaled  values  are  typically  stored  back  in  bins  <  n  <  N 
and  plotted  at  the  upper  end  of  the  interval  (0,2ii//\^).  However, 
this  upper  end  of  the  plot  must  be  recognized  as  the  estimate  of 
EDF  E(u)  just  to  the  left  of  the  origin,  that  is,  u  <  0. 

Usually,  only  the  values  for  0  <  n  <  are  of  interest  because 
they  encompass  the  small  tail  probability  values  of  E(u). 

The  infinite  sum  on  m  in  equation  (24)  must  be  truncated 
where  the  magnitude  of  complex  ratio  “  !!■)/(  r  + 

becomes  insignificant.  Values  of  this  complex  ratio  for  m  >  N 
are  simply  additively  prealiased  into  bin  m  MODULO  N,  thereby 
avoiding  truncation  error.  Notice  that  FFT  size  N  has  no  effect 
on  accuracy;  it  merely  controls  the  spacing  =  2ii/(NA^)  of 
output  values  in  equations  (24)  and  (25). 


18 


EXAMPLES 


Four  different  characteristic  functions  are  considered  in 
this  section.  For  each  CF  f(^),  the  corresponding  PDF  p(u)  is 
calculated  and  plotted,  along  with  an  error  estimate.  In  the 
second  half  of  this  section,  the  corresponding  EDFs  E(u)  are 
evaluated,  in  addition  to  an  error  estimate.  Deeper  tail 
probabilities  than  considered  in  these  examples  are  possible 
through  larger  choices  of  the  shift  parameter  r.  All  the  FFT 
sizes  were  taken  at  N  =  1024  in  this  section. 


CHI-SQUARE  PDF 

The  CF  of  interest  is  f{^)  =  (1  -  which  corresponds 

to  a  random  variable  x  that  is  composed  of  a  sum  of  20  indepen¬ 
dent,  squared  Gaussian  random  variables  with  zero  mean  and  common 
variance  Because  x  >  0,  additive  constant  c  is  taken  as  zero. 

Consider,  first,  the  standard  FFT  approach  (1)  on  this  CF, 
that  is,  with  shift  parameter  r  =  0.  The  resultant  PDF  p(u)  is 
displayed  in  figure  1,  where  sampling  increment  =  2jt/80  was 
used  in  equation  (7)  and  sequel.  The  PDF  p(u)  decays  until  it 
reaches  the  inherent  round-off  noise  of  the  computer,  which  is 
roughly  E-15  relative  error  (64-bit  representations).  Lower 
probability  (density)  estimates  are  not  possible  for  this  example 
with  this  standard  FFT  approach.  The  maximum  abscissa  (period) 
on  this  plot  is  2ii/A^,  which  is  80  for  the  above  choice  of  A^. 
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Figure  1.  Chi-Square  p(u)  via  Standard  FFT  Technique 


Observe  that  p(u)  drops  from  E-6  to  E-16  in  a  span  of  values 
from  u  =  32  to  u  ”  61,  that  is,  a  span  of  29  in  u.  To  reach 
probability  (density)  level  E-76,  for  example,  an  additional  span 
of  approximately  6  *  29  =  174  in  u  is  required.  To  avoid 
aliasing,  a  new  sampling  increment  satisfying  61  +  174  =  2n/A^, 
namely,  =  2ii/235,  is  required.  That  is,  a  noise-free  computer 
would  have  required  A^  =  2ji/235  to  get  to  level  E-76  by  the 
standard  technique  with  r  =  0. 
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On  the  other  hand,  using  shift  parameter  value  r  =  0.7,  along 
with  c  =  0  and  =  2ii/240,  the  periodic  function  p(u),  defined 
in  equation  (7),  is  presented  in  figure  2.  Naturally,  p(u)  has  a 
relative  error  about  E-15;  furthermore,  this  small  sampling 
increment  is  observed  to  be  necessary  to  avoid  aliasing  at  the 
right  tail  of  p(u). 

When  the  results  for  p(u)  are  used  in  equation  (16),  the 
resultant  estimate  for  the  desired  PDF  p(u)  is  as  shown  in  figure 
3.  The  superposed  dashed  line  is  an  error  estimate,  starting 
from  E-12  near  the  top  left  and  going  to  the  noisy  region  at  the 
bottom  right  end.  The  relative  error  of  the  estimated  p(u) 
varies  over  the  range  of  u  values,  gradually  deteriorating  to 
useless  values  for  u  >  200.  However,  approximately  six 
significant  decimal  digits  are  still  available  at  the  E— 50 
probability  level.  For  this  example,  the  exact  PDF  p(u)  is 
available  and  has  also  been  plotted  in  figure  3;  it  overlies  the 
estimated  PDF  until  the  probability  values  drop  below  E-75. 

GAUSSIAN  PDF 

2 

The  normalized  Gaussian  CF  is  f(^)  =  exp(-^  /2  +  ic^),  where 
additive  constant  c  is  taken  as  2.  With  shift  r  =  7  and  sampling 
increment  A^  =  2ii/20,  the  periodic  function  p(u)  is  as  displayed 
in  figure  4.  The  sampling  increment  could  not  have  been  taken 
much  smaller  without  incurring  aliasing.  Also,  the  constant  c 
is  at  its  minimum  value,  just  keeping  positive  values  at  u  =  0+. 


21 


u 


Figure  3. 


Chi-Square  p(u) 


The  corresponding  estimated  PDF  p(u)  is  plotted  in  figure  5, 
along  with  the  estimated  error.  The  estimated  error  line  joins 
the  two  noisy  regions  at  each  end  of  the  interval  (0,20),  which 
are  easily  located  in  p(u)  in  figure  4.  For  this  example,  there 
are  approximately  seven  decimals  of  accuracy  at  the  E-40 
probability  level.  The  exact  PDF  p(u)  overlies  the  estimated  PDF 
except  at  the  edges  of  the  interval  (0,20);  this  type  of  result 
is  again  anticipated  by  the  plot  of  figure  4. 


BRANCH  POINT  PDF 

The  CF  for  this  example  has  multiple  branch  points,  namely, 

10  -k 

f(^)  =  rT(i  - 

n=i 

This  CF  corresponds  to  the  sum  of  10  independent,  squared  zero- 
mean  Gaussian  random  variables  with  different  variances,  l/(2n). 
The  corresponding  PDF  p(u)  is  not  available  in  closed  form. 

The  relatively  slow  decay  of  this  CF,  namely,  as  1900/^^, 
forces  numerous  evaluations  of  equation  (27)  to  be  undertaken. 

In  fact,  with  increment  =  2ii/200,  over  160,000  complex  samples 
were  required  until  |f( ^) |  became  sufficiently  small  to  control 
the  truncation  error.  Nevertheless,  employment  of  prealiasing 
operation  (12)  enabled  storage  and  execution  of  only  an  N  =  1024 
point  FFT  according  to  equation  (13). 
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Figure  4.  Gaussian  p(u) 


The  function  p(u)  resulting  from  operation  (13)  is  displayed 
in  figure  6  for  shift  r  =  0.8.  The  corresponding  estimated  PDF 
p(u)  is  given  in  figure  7,  along  with  the  estimated  error.  The 
result  for  p(u)  has  approximately  six  digits  of  significance  at 
the  E-40  probability  level. 


ESSENTIAL  SINGULARITY  PDF 

The  last  characteristic  function  of  interest  is 

,10  1  N  10 

f(^)  =  expfi^  C  i 

''  n=l  n=l 

This  CF  has  branch  points  and  essential  singularities  at  ^  =  -in 
and  corresponds  to  a  sum  of  10  independent,  squared  Gaussian 
random  variables  with  a  common  mean  of  1  and  different  variances 
l/(2n). 

The  exponential  in  equation  (28)  causes  a  more  rapid  decay 
with  in  fact,  only  200  samples  of  the  CF  were  necessary  at 
sampling  increment  =  2n/200,  using  shift  r  =  0.6,  to  control 
the  truncation  error.  The  resultant  p(u),  plotted  in  figure  8, 
again  has  relative  error  in  the  E-15  range.  The  corresponding 
estimated  PDF  p(u)  and  its  estimated  error  are  displayed  in 
figure  9.  There  are  approximately  five  digits  of  significance  at 
the  E“40  probability  level. 
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CHI-SQUARE  EOF 


—  10 

The  CF  of  interest  is  again  f(^)  =  (1  -  i^)  •  However,  for 

the  following  examples,  the  quantity  of  interest  is  the  EOF  E{u). 
For  sampling  increment  =  2ii/240,  shift  r  ==  0.6,  and  additive 
constant  c  =  0,  the  periodic  function  E(u)  of  equations  (21)  and 

(24)  is  presented  in  figure  10.  The  nonzero  values  of  E(u)  just 
below  u  =  240  are  the  nonzero  values  of  E(u)  for  u  <  0  being 
periodically  repeated;  this  exponential  decay  of  E(u),  according 
to  equation  (20),  shows  up  as  a  straight  line  on  the  logarithmic 
ordinate.  The  relative  error  of  E(u)  in  figure  10  in  the  noisy 
region  is  E-15. 

For  this  example,  take  the  noise-region  integer  of  equation 

(25)  as  Nj  =  0.75  N  =  768.  Then,  the  estimated  EDF  E(u), 
obtained  from  equations  (25)  and  (26),  is  presented  in  figure  11. 
For  u  less  than  150,  the  estimated  and  exact  EDFs  overlap.  Also, 
for  200  <  u  <  240,  the  estimated  E(u)  recovers  its  correct  value 
of  1;  recall  that  this  region  corresponds  to  u  <  0  in  the 
original  E(u) . 

The  accuracy  of  estimated  EDF  E(u)  is  about  five  decimals  at 
the  E-40  probability  level.  The  error  estimate  (dashed  line)  is 
obtained  by  drawing  a  straight  line  between  the  noisy  regions 
established  in  figure  10;  the  white-noise  level  in  figure  10  is 
translated  by  equations  (23)  or  (25)  into  a  line  with  slope 
-r  log(e)  =  -0.434  r  on  the  logarithmic  ordinate  of  figure  11. 
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Figure  10.  Chi-Square  E(u) 


Figure  11-  Chi-Square  E(u) 


GAUSSIAN  EDF 


2 

The  CF  for  this  example  is  f(^)  =  exp(-^  /2  +  ic£;)  ,  where 
constant  c  is  taken  as  6.  For  increment  *=  2n/24  and  shift 
r  -  5,  the  periodic  function  E(u)  is  as  displayed  in  figure  12. 
For  noise-region  integer  =  7/8  N  =  896,  the  corresponding 
estimated  E(u)  is  plotted  in  figure  13.  The  stability  is 
approximately  six  decimals  at  the  E— 30  probability  level.  The 
estimated  E(u)  never  recovers  its  unity  value  for  u  just  less 
than  24  because  the  periodic  function  in  figure  12  was  not 
noise-free  in  that  region  (in  contrast  to  figures  10  and  11). 


BRANCH  POINT  EDF 

The  CF  f(^)  was  given  in  equation  (27).  The  increment 
must  be  taken  smaller  for  the  EDF  than  for  the  PDF  because  of  the 
"left"  tail  of  E(u)  for  u  <  0;  see  equation  (20).  On  the  other 
hand,  the  number  of  samples  decreased  to  about  60,000  because  of 
the  extra  1/K  decay  factor  in  equations  (18)  and  (21).  The 
periodic  function  E(u)  is  given  in  figure  14.  The  left  tail  of 
E(u)  shows  up  in  the  region  200  <  u  <  240. 

The  corresponding  estimated  EDF  E(u)  is  displayed  in  figure 
15.  It  has  seven  significant  digits  at  the  E-40  probability 
level.  The  error  estimate  has  slope  -0.434  r,  which  is  -0.35  for 
this  example  where  r  =  0.8  has  been  used. 
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Figure  15.  Branch  Point  E(u) 
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ESSENTIAL  SINGULARITY  EOF 


The  CF  for  this  example  is  given  by  equation  (28).  The 
exponential  decay  of  f(?)  allows  accurate  calculation  of  the 
periodic  function  E(u)  with  just  over  200  samples  of  the  CF .  The 
left  tail  in  figure  16,  just  below  u  =  240,  is  again  the  periodic 
repetition  of  the  exponential  behavior  of  E(u)  for  u  <  0,  as 
indicated  in  equation  (20). 

The  corresponding  estimated  EOF  E(u)  is  given  in  figure  17, 
along  with  an  error  estimate  that  has  a  slope  of  -0.434  *  0.6  = 
-0.26  for  r  =  0.6  in  this  example.  There  are  approximately  four 
digits  of  significance  in  E(u)  at  the  E-40  probability  level. 
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Figure  16.  Essential  Singularity  E(u) 


Figure  17.  Essential  Singularity  E(u) 
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SUMMARY 


An  efficient,  fast,  and  accurate  FFT  technique  for  obtaining 
small  tail  probabilities  for  both  the  PDF  and  the  EDF ,  directly 
from  the  CF,  has  been  derived  and  demonstrated  numerically  for 
several  examples.  The  method  is  especially  useful  when  analytic 
or  asymptotic  expressions  for  the  probabilities  are  unavailable 
or  unknown. 

By  choosing  the  shift  parameter  r  closer  to  the  highest 
singularity  of  CF  f(^)  in  the  complex  ^-plane,  smaller  values  of 
the  tail  probabilities  of  the  PDF  and  EDF  can  be  realized.  The 
cost  in  this  approach  is  that  the  sampling  increment  must  then 
be  taken  smaller  to  avoid  aliasing.  Although  taking  a  smaller 
necessitates  more  computer  time  and  effort,  it  does  not  have  to 
involve  more  storage  because  prealiasing  can  be  advantageously 
employed  to  keep  the  FFT  size  N  at  very  reasonable  values.  The 
FFT  size  N  has  no  effect  upon  the  errors  caused  by  aliasing  and 
truncation;  rather,  N  merely  controls  the  spacing  at  which  the 
PDF  and  EDF  outputs  are  calculated. 

A  sample  BASIC  program  for  calculating  the  PDF  p(u)  is  listed 
in  appendix  A,  and  the  corresponding  program  for  the  EDF  E(u)  is 
given  in  appendix  B.  Inputs  required  of  the  user  are  additive 
constant  c,  shift  r,  sampling  increment  A^,  and  the  particular  CF 
f<^).  Some  trial  and  error  may  be  necessary  in  selecting  the 
best  values  of  r  and  A^. 
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APPENDIX  A  —  BASIC  PROGRAM  FOR  PDF 


10 

SMALL  TAIL  PROBABILITIES 

"CHI-SQUARE-PDF" 

20 

R  = 

.7  ! 

AMOUNT  OF  CONTOUR  SHIFT  USED 

30 

Dx 

=2.*PI/240.  ! 

SAMPLING  INCREMENT  IN  X 

40 

C  = 

0.  ! 

ADDITIVE  CONSTANT 

50 

N= 

1024  1 

SIZE  OF  FFT 

60 

Y1 

=  -90.  ! 

LOWER  LIMIT  OF  LGT  RANGE 

70 

Y2 

=  0.  ! 

UPPER  LIMIT  OF  LGT  RANGE 

80 

To1=l.E-36  ! 

TOLERANCE  ON  NAG  SQ  CHAR  FN 

90 

N1 

=N-1 

100 

REDIM  Co£<0: N/4) , X<0: N1 > , Y<0: HI ) 

110 

DIM  Cos- <4096)  ,  X<  16384)  ,  Y<  16384) 

120 

DOUBLE  H,Nl,K,M,Ns  ! 

INTEGERS,  NOT  DOUBLE  PRECISION 

130 

fl= 

2.*PI/N 

140 

FOR  K=0  TO  N/4 

150 

Co 

s<K)=C0S<A#K)  ! 

QUARTER-COSINE  TABLE  IN  Co=<«) 

160 

NEXT  K 

170 

fil=10.  i 

CF  f<xi)  =  1/<1  -  i  >:i)-''alpha 

180 

Rl=l.-R 

190 

F0=1./R1‘^R1 

200 

X<0>=.5 

210 

Y<0>=0. 

220 

M  =  0 

230 

M  =  M+1 

240 

X=Dx*M 

250 

CfiLL  Poujer<Rl,-X,-fll,fir,fli> 

260 

Fr=fir/F0  ! 

REAL  CF 

270 

Fi=fii/F0  1 

IMAG  CF 

280 

K=M  MODULO  N  ! 

PREALIASING 

290 

X<:K>  =  XCK)+Fr 

300 

Y<K)=YCK)+Fi 

310 

IF  (Fr*Fr+Fi *Fi >Tol >  THEN  230 

320 

CALL  Ff t 14<N, Co£<*) , X<*) , Y<*) ) 

-"'330 

MAT  X=X»<Dx/PI) 

340 

Big=MRX<X<*)) 

350 

GINIT 

360 

PLOTTER  IS  "GRAPHICS" 

370 

GRAPHICS  ON 

330 

WINDOW  0,N,-18,0 

390 

LINE  TYPE  3 

400 

GRID  N/8,2 

410 

LINE  TYPE  1 

420 

FOR  Ns=0  TO  N1 

430 

Pn=X<Ns)  ! 

PDF  TILDA 

440 

IF  Pn<>0.  THEN  470 

450 

PENUP 

460 

GOTO  480 

470 

PLOT  Ns,LGT<RBS<Pn)/Bi g)  ! 

NORMALIZED  PDF  TILDA 

480 

NEXT  Ns 

490 

PENUP 

500 

PAUSE 

A-1 


!  U  INCREMENT 


510  Du=2.  *PI/<N*Dx) 

520  FOR  Ns=0  TO  N1 

530  U=Du»Ns 

540  X<:Ns)=F0*EXP<-R*IJ)*X'::Ns.>  !  PDF  ESTIMATE 

550  NEXT  Ns 

560  GCLEAR 

570  WINDOW  0,N,Y1,Y2 

580  LINE  TYPE  3 

590  GRID  N/8, 10 

600  fil=fil-l. 

610  Ga=FNGamma< A 1 ) 

620  FOR  Ns=l  TO  N  STEP  10 

630  U=Du*Ns 

640  Pe=FNEx<U-Al*LOG(U) )/Ga  !  PDF  EXACT 

650  PLOT  NSjLGTCPe) 

660  NEXT  Ns 

670  PENUP 

680  PAUSE 

690  LINE  TYPE  1 

700  FOR  Ns=0  TO  N1 

710  Pn=X<Ns)  !  PDF  ESTIMATE 

720  IF  Pn<>0.  THEN  750 

730  PENUP 

740  GOTO  760 

750  PLOT  Ns,LGT<RBS<Pn)> 

760  NEXT  Ns 

770  LINE  TYPE  3 

780  MOVE  0,LGT<AES<X<0)))  !  ERROR 

790  DRAW  N1 , LGTCABS<X<N1 ) ) )  !  ESTIMATE 

800  PENUP 

810  PAUSE 

820  END 

830  ! 

840  DEF  FNEx<X)  !  EXP<-X) 

850  IF  X>708.  THEN  RETURN  0. 

860  RETURN  EXP(-X) 

870  FNEND 

880  ! 

890  SUE  PouerCX,  Y.Real  ,U,V>  !  PRINCIPAL  POWER  Z'-Real 

900  T=X*X+Y*Y 

910  IF  T>0.  THEN  950 

920  -  U=V=0. 

930  IF  Real=0.  THEN  U=l. 

940  SUBEXIT 

950  F=EXP<.5*Real*L0G<T)) 

960  IF  X=0.  THEN  A=. 5*PI*SGN<Y) 

970  IF  XO0.  THEN  A  =  ATN<Y/X> 

980  IF  X<0.  THEN  A  =  A  +  PI*(:i . -2.  Y<0.  )  ) 

990  T=Real*A 

1000  U=F*COS<T> 

1010  V=F*SIN(T> 

1020  SUBEND 

1030  ! 

1040  DEF  FNGanuftat X )  !  HART,  paQe-  135,  #5243! 

1370  ! 

1380  SUE  Fftl4CIiOUBLE  N,REAL  Cos < *  > ,  X ( * >  ,  Y ( * >  )  !  H<  =2' 1  4=  1  6334  ;  0  SUES 
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APPENDIX  B  -  BASIC  PROGRAM  FOR  EDF 


10  I 
20 
30 
40 
50 
60 
70 
80 
90 
100 
110 
120 
130 
140 
150 
160 
170 
180 
190 
200 
210 
220 
230 
240 
250 
260 
270 
280 
290 
300 
310 
320 
330 
340 
350 
360 
370 
380 
390 
400 
410 
420 
430 
440 
450 
460 
470 
480 
490 
500 


SMfiLL  TfilL  PROBRBILITIES  " CH I -SQURRE-EnF " 

Dx  =  2.*PI.^240.  !  SRriPLIHG  INCREMENT  IN  X 

R=.6  !  RMOUNT  OF  CONTOUR  SHIFT  USED 

N=1024  !  SIZE  OF  FFT 

Yl=-70.  !  LONER  LIMIT  OF  LGT  RRNGE 

Y2=10.  !  UPPER  LIMIT  OF  LGT  RANGE 

Tol=l.E-30  !  TOLERRNCE  ON  MAG  SQ  CHAR  FN 

N1=N-1 

REDIM  Cos < 0 : N/4 > , X < 0 : N 1 ) , Y < 0 :  N 1 ) 

DIM  Cos<4096) , X( 16384) , Y( 16384) 

DOUBLE  N,Nl,K,M,Ns,Nf  !  INTEGERS,  NOT  DOUBLE  PRECISION 

fi=2. *PI/N 

FOR  K=0  TO  N/4 

Cos<K)=COS(A*K)  !  QUARTER-COSINE  TABLE  IN  r;ns<*) 

NEXT  K 

MAT  X=<0. ) 

MAT  Y=<0. ) 

81=10-  !  CFfCxi)  =  l/<l-i  xi)"*alph5L 

R1=1.“R 

R2  =  R*R 

X<0)=. 5/<Rl^Rl *R) 

M  =  0 
M  =  M+1 
X  =  Dx*M 

CALL  PowerCRl,-X,-Al , Ar,Ri ) 
fi=R2+X*X 

Fr=<flr*R+Ai *X)/fi  !  REAL  CF 

Fi =<Ri *R“Rr*X)/fi  !  IMAG  CF 

K=M  MODULO  N  I  PRERLIRSING 

X<K)=X<K)+Fr 
Y<K)=Y<K)+Fi 

IF  <Fr*Fr+Fi*F1 ) >Tol  THEN  230 
CALL  Fftl4<N,Cos<«'),X<*),  Y<*)) 

MAT  X=X*<Dx/PI) 

Big=MRXCX(*)) 

GINIT 

PLOTTER  IS  “GRAPHICS” 

GRAPHICS  ON 
NINDON  0,N,-18,0 
LINE  TYPE  3 
GRID  N/8,2 
LINE  TYPE  1 
FOR  Ns=0  TO  N1 

En=X<Ns)  !  EDF  TILDA 

IF  En<>0.  THEN  480 
PENUP 
GOTO  490 

PLOT  Ns,  LGT<RBS<Eri)/Bi  g)  !  NORMALIZED  EDF  TILDA 

NEXT  Ns 

PENUP 


510 

520 

530 

540 

550 

560 

570 

580 

590 

600 

610 

620 

630 

640 

650 

660 

670 

680 

690 

700 

710 

720 

730 

740 

750 

760 

770 

780 

790 

800 

810 

820 

^830 

"840 

850 

860 

870 

880 

890 

900 

910 

920 

930 

940 

950 

960 

1100 

1110 

1440 

1450 


INPUT  "INPUT  FRPCTION  OF  PERIOD", Frac 
Nf=N*Frac  ,  , 

.Du=2.*PI/<N*Dx)  !  U  INCREMENT 

FOR  Ns=0  TO  Nf 
U=Du^Ns 

Xc:Ns)=EXP(-R*U)^X<Ns)  !  EDF  ECu) 

NEXT  Ns 
Rp  =  R*2.  '^PI/Dx 
FOR  Ns=Nf+l  TO  N1 
U=Du*Ns 

XCNs)=EXP<-R4^U  +  Rp)*X<Ns)  !  EDF  E<u) 

NEXT  Ns 
GCLERR 

NINDOW  0,N,Y1,Y2 
LINE  TYPE  3 
GRID  N/8, 10 
fll=Rl-l. 

Ga=FNGamnia(fll  > 

FOR  Ns=0  TO  N  STEP  10 

U  =  Du*Ns 

S=T=1. 

FOR  K=1  TO  Fll 
T=T#U/K 
S  =  S  +  T 
NEXT  K 

Ee=EXP(-U>^S  !  EDF  EXRCT 

IF  Ee>0.  THEN  800 

PENUP 

GOTO  810 

PLOT  Ns,LGT<Ee) 

NEXT  Ns 
PENUP 

LINE  TYPE  1 
FOR  Ns=0  TO  N1 

En=X(Ns>  !  ESTIMRTED  EDF 

IF  En<>0.  THEN  890 

PENUP 

'  GOTO  900 

PLOT  Ns,LGT<RES(En)> 

NEXT  Ns 

PENUP 

PfiUSE 

GRRPHICS  OFF 
END 

I 

SUE  Power(X, Y,Rea1 ,U, V)  !  PRINCIPRL  POWER  Z-Real 
! 

DEF  FNGamrfiac:X)  !  HRRT,  page  135,  #5243 

j 

SUE  Fftl4<D0UELE  N,RERL  Cos  <  *  )  ,  X < *  >  ,  Y  <  *  )  )  !  H<  =  2'^- 1  4  =  1  6384  ;  0  SUES 
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